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Abstract —Because of the attractiveness of the canonical 
polyadic (CP) tensor decomposition in various applications, 
several algorithms have been designed to compute it, but efficient 
ones are still lacking. Iterative deflation algorithms based on 
successive rank-1 approximations can be used to perform this 
task, since the latter are rather easy to compute. We first present 
an algebraic rank-1 approximation method that performs better 
than the standard higher-order singular value decomposition 
(HOSVD) for three-way tensors. Second, we propose a new iter¬ 
ative rank-1 approximation algorithm that improves any other 
rank-1 approximation method. Third, we describe a probabilistic 
framework allowing to study the convergence of deflation CP 
decomposition (DCPD) algorithms based on successive rank-1 
approximations. A set of computer experiments then validates 
theoretical results and demonstrates the efficiency of DCPD 
algorithms compared to other ones. 

Index Terms —rank-1 approximation, Canonical Polyadic, ten¬ 
sor decomposition, iterative deflation, blind source separation 

I. Introduction 

In the last years, tensors have been playing an important role 
in many applications such as blind source separation fl|, lf2), 
telecommunications 11, chemometrics 0, neuroscience (3), 
sensor array processing (6) and data mining 0. The attractive¬ 
ness behind tensors lies in the uniqueness of their canonical 
polyadic (CP) decomposition under mild conditions (8), which 
is a powerful property not shared by standard matrix-based 
tools. There are several methods to compute the CP tensor 
decomposition. We will point out here some of the most used 
methods among many others. For the exact CP decomposition, 
0 proposes a direct computation to decompose 2 x n x n 
tensors. In flOfl . a generalization of Sylvester’s algorithm is 
described for decomposing symmetric tensors. In ifTTl . one 
can use simultaneous matrix diagonalization by congruence, 
provided that the rank of the tensor is smaller than its greatest 
dimension. An approach based on eigenvectors of tensors is 
proposed in fl2l . 

In practice, tensors are corrupted by noise so that one needs 
to compute an approximate decomposition of given rank. 
Computing the exact CP decomposition is difficult fl3l , but 
finding a lower-rank approximation is even harder. In fact, 
this is an ill-posed problem in general m. Nevertheless, 
some useful algorithms have been conceived to solve locally 
the low-rank approximation problem. This kind of algorithms 
can be found in m, m, im, ei, Ei- One of the most 
widely used is the alternating least squares (ALS) algorithm 
0, which is an iterative method that consists in conditionally 
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updating in an alternate way, the matrix factors composing the 
CP decomposition. Other gradient and Newton-based methods 
estimate the factor matrices all-at-once. Howsoever, all theses 
algorithms have disappointing convergence properties O, 
[201]. Another kind of algorithms is based on rank-1 deflation. 
It is known that the conventional deflation works for matrices 
but not for tensors 121) . In | |22| , the authors propose deflation 
methods that work only if the rank of the tensor is not larger 
than its dimensions. In l23l . ALS is used to update the columns 
of matrix factors in a deflation procedure, but for non-negative 
tensors. However, these deflation methods strongly depend on 
initialization and no convergence study has been conducted. 

In the same vein of iterative deflation, the authors pro¬ 
posed in lf24l a deflation-based CP decomposition (DCPD), 
based on successive rank-1 approximations computed by low- 
complexity methods. Other rank-1 approximation methods can 
be used in DCPD, for instance, ALS. However, the latter 
exhibits an unbounded complexity and no satisfactory conver¬ 
gence study is available for rank-1 approximation apart l26fl . 
which shows results on the global convergence for generic 
tensors in the sense that, for any initialization, ALS con¬ 
verges to the same point in general. A quasi-Newton method 
defined on Grassmannians are developed in El- However, 
they exhibit an unbounded complexity to compute rank-1 
approximations, since these methods need to be iterated. In 
[28), the best rank-1 approximation problem can be computed 
by means of an algebraic geometry moment method, but 
it is only applicable to very small tensor dimensions since 
the number of variables grows exponentially when building 
convex relaxations. Moreover, even for small dimensions, its 
convergence is very slow. In (29) . the authors propose an 
improvement method for |28| based on border basis relaxation, 
but again the method is limited to small tensor dimensions. 
Semidefinite relaxations are proposed in 130) to compute rank- 
1 approximations, however the convergence becomes very 
slow when dimensions are large. 

In this paper, we report mainly three contributions. First, 
we propose an algebraic rank-1 approximation method, namely 
the sequential rank-1 approximation and projection (SeROAP), 
which can perform better than the standard truncated high- 
order singular value decomposition (THOSVD) (25). Indeed, 
we prove that the rank-1 approximation performed by SeROAP 
is always better than the one obtained from THSOVD for 
three-way tensors. Moreover, for large dimensions and small 
orders, we show that the computational complexity of SeROAP 
is dramatically lower than that of THOSVD. 

Second, we propose an alternating eigenvalue rank-1 iter¬ 
ative algorithm for three-way tensors, namely CE (coupled- 
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eigenvalue), that improves other rank-1 approximation algo¬ 
rithms. We prove that if the solution obtained from some rank- 
1 approximation algorithm (e.g., SeROAP, THOSVD) is the 
input of the CE algorithm, the performed rank-1 approximation 
remains the same in the worst case. We also prove that 
the convergence to a stationary point is always guaranteed. 
Actually, results have shown that when the initialization of 
the CE algorithm is close enough to the global solution, it 
recovers the best rank-1 approximation. Furthermore, when 
one dimension is much larger than the other two dimensions, 
the computational complexity of the CE algorithm can be 
lower than that of the standard ALS algorithm. 

Third, we perform a theoretical study on deflation in order 
to analyze the convergence of the DCPD algorithm. In a first 
stage, we show that the norm of residuals is monotonically 
reduced within the iterative deflation process. We also prove 
that the DCPD algorithm recovers the exact CP decomposition 
of a given tensor when residuals do not fall within a cone with 
an arbitrary small volume. In a second stage, we prove that the 
iterative deflation method can reduce the norm of the initial 
residual by a factor smaller than (sin(/3)) L ~ 1 (/? being the 
angle of a suitable cone where the residuals can fall in) after 
L iterations with high probability, when tensors are distributed 
according to an absolutely continuous probability measure, and 
the probability function of residuals is continuous on some 
suitable angular interval. We also present a conjecture stating 
the existence of probability measures ensuring the convergence 
of the DCPD algorithm to an exact decomposition with high 
probability. 

The paper is organized as follows. In Section [TIT] some 
standard iterative methods and the DCPD algorithm as well as 
SeROAP and THOSVD methods are described. The computa¬ 
tional complexity per iteration for each algorithm is provided. 
The next two sections form the core of the paper. In Section 
IIV1 we first prove that SeROAP performs better than THOSVD 
as far as rank-1 approximations of three-way tensors is con¬ 
cerned. Then in a second part, the CE algorithm is presented 
as well as the proof that it can refine any other rank-1 approxi¬ 
mation method and the proof of its convergence. In Section [VI 
among other theoretical results, we study conditions ensuring 
the convergence of the DCPD algorithm. Finally, in Section 
ED computer results show satisfactory performances of the 
proposed DCPD and rank-1 approximation methods, compared 
to other related algorithms, even under noisy scenarios. 

II. Notation 

The notation employed in this paper is the following. Scalar 
numbers are denoted by lowercase letters and vectors by 
boldface lowercase ones. For matrices and tensors, boldface 
capital and calligraphic letters are used, respectively. Plain 
capitals are used to denote array dimensions. Depending on 
the context, greek letters can denote either scalars, vectors, 
matrices or tensors. The symbols □, 0, S3 and 0 denote 
the Hadamard, Khatri-Rao, Kronecker and tensorial products, 
respectively, and + denotes matrix pseudo-inversion. The Eu¬ 
clidean scalar product between tensors is denoted by (X. Li) = 
T n-iN U t 1 ...i N - The angle between two tensors U and 


V will refer to arccos{|(ZX V)|}/||W||f ||V||f G [0, 7t/2]. 
|| • || then denotes the Frobenius norm induced by the previous 
scalar product. We shall also use the l 2 operator norm || • ||2 
for matrices, which corresponds to the largest singular value. 
The mode-n unfolding of a tensor X is denoted as T( n y as 
proposed in fl6l . 7~(i,:,k) and denote 

vector slices of tensor X. 

The operator vec is the vectorization operator that stacks the 
columns of a matrix into a long column vector, and Unvec 
is its reverse operator. C k is the set of functions having 
continuous /cth derivatives. Finally, IK is either the real or the 
complex field. 

III. Description of Algorithms and Complexity 
Analysis 

This section presents the description of some CP decom¬ 
position algorithms. In order to support further results, the 
complexity per iteration is calculated here for each algorithm 
using Landau’s notation, denoted by O , and counting only 
multiplications, as recommended in 1311 . From Section llll-AI 
to Section 1III-B1 we describe two classical algorithms known 
in literature: ALS and conjugate gradient ns. In Section 
IIII-CI the DCPD algorithm is presented. 

For the CP decomposition algorithms described in the 
following, the input parameter R denotes the rank of the output 
tensor. Assuming Rq is the rank of the input tensor X, if 
R 0 < R , then the algorithms perform an exact decomposition. 
On the other hand, if Rq > R, a lower rank-f? approximation 
is computed. 

A. Alternating least squares (ALS) 

The most commonly used algorithm for solving the CP 
decomposition is ALS S). The goal is to update alternately 
each factor matrix in each iteration by solving a least squares 
problem conditioned on previous updates of the other factor 
matrices. There is no guarantee of convergence to the global 
solution, nor even to a critical point. The implementation is 
quite simple and it is detailed in Alg{T] 


input : X € K 7lX/2X "' x/,v : input data, 

R: rank parameter. 

output: A (n ) e K InXR , for n = 1,... ,TV: factor 
matrices 

Initialize A (1 ^,..., A (jV) 

repeat 

for n = 1 to N do 

V «- A il)J A {1) □ • • ■ □ A^-^ □ 

A (n+ l)T A (n+l) H H A (N)T A (N) 

A (n) = T (n) (A (iV) © • • • © A (ra+1) © A (n " 1} © 
•■•© A {1) )V + 

end 

until some stopping criterion is satisfied ; 

Algorithm 1: ALS algorithm 


The complexity per iteration (repeat loop) of ALS may be 
calculated as follows. The computation of matrix V needs 
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(Ji + ... + I n -1 + In+i ■ ■ ■ In)R 2 + (N - 2)R 2 operations 
(multiplications). The pseudo-inverse of V is calculated by 
resorting to an SVD. For a roxn rank-r matrix, with m > n > 
r, the explicit calculation of diagonal, left singular, and right 
singular matrices require 2 mn 2 — 2n 3 /3, 5mr 2 — r 3 /3, and 
5nr 2 —r 3 /3 multiplications, respectively f3Tl . Thus, assuming 
for simplicity that V is a non-singular matrix, the number of 
operations to calculate its pseudo inverse is 35f? 3 /3 + R 2 . For 
updating A 1 n \ we need 

N N 

R li + R Ij + RR 2 

j=l 

multiplications. These calculations must be performed for each 
n £ {1,..., N}. Thus, the number of operations per iteration 
for ALS is dominated by the term composed of the product 
of all dimensions. Hence, 


#op 



B. Conjugate gradient (CG) 

The conjugate gradient algorithm (CG) is a faster algorithm 
than the well known gradient descent im Here, we use the 
optimum step size and the Polak-Ribiere implementation l32l 
for updating the parameter /3 in the algorithm presented in 
AlgE] The number of operations for computing each vector 
9n is 


N N N 

R 2 r 3 +(N- 2 )R 2 + InR 2 + R n I ’- + R n L i 

j —1 J —1 

The computation of the step size //* is dominated by the 
number of multiplications needed to determine all the coef¬ 
ficients but one of a 2A’-degree polynomial generated from 
the enhanced line search (ELS) method, which is given by 

03 : 

N 

(2 N R + 0{N 2 })l[l J 

j =i 

The parameter ft requires only 1 + 2RYlf=i^j multi¬ 
plications. Hence, the CG algorithm with ELS has a total 
complexity given by 


#op = O j ((2 


N 


N 


1)R + N 2 )1[Ij 


i=i 


( 2 ) 


input T € K IlXl2X "' xlN : input data, 

R: rank parameter. 

output: A (n) € K / " x - r , for n = 1,..., TV: factor 
matrices 

Initialize A (1) ,..., A (Ar T 
p x— [uec T (A (1 ^) • • • uec T (A (Ar ' ) )] T ; 

repeat 

for n = 1 to N do 

Compute the gradient g, n ^ with respect 

uec(A (n) ); 

end 

a <- [g T {1) ■ ■ ■ gJ N) } J -, 

if n = 1 then 

I d i - g\ 

end 

Compute the optimum step size 
Update /3 according to Polak-Ribiere; 

A p x— /i*d; 

p<-p + A p; 

d i — —g [Id\ 

until some stopping criterion is satisfied ; 

Extract A"-* from p , for n = 1,..., N. 

Algorithm 2: CG algorithm 


1) Truncation of higher order singular value decomposi¬ 
tion - THOSVD: The algorithm is described in A IgU] For 
computing the first right singular vector, we do not need 
compute the complete SVD. According to lf33l . we can 
compute the best rank-1 approximation of a m x n matrix 
in k steps by using the Lanczos algorithm, with a complexity 
0{2kmn}. Hence, the accumulated complexity computed for 
all u n is equal to 0{2NkY[^ = i Ij}- The computation of li 
requires IljLi Ij fl°P s - The contraction to obtain A also needs 
Ij operations. To sum up, the total number of operations 
of THOSVD is of order: 

r N 

ol ( 2 Nk + 2 ) Y[ Ij 

j=i 


input : T G K IlXl2X '" xlN : input data 

output: X e K /lX/2X - X A : rank-1 approximation 

for n = 1 to N do 

] u n x— first left singular vector of T( n ) ; 

end 

U t- ®n=l U n\ 

A <- <T,M); 

X ^ X U. 

Algorithm 3: THOSVD algorithm 


C. Deflation-based CP decomposition (DCPD) 

The computation of a rank-1 approximation is the key of the 
DCPD algorithm. We present here two methods for computing 
a rank-1 approximation referred to as THOSVD and SeROAP 

ED- 


2) Sequential rank-1 approximation and projection - 
SeROAP: Without loss of generality, consider I\ > Ij > ■ ■ ■ > 
In- The SeROAP algorithm ll24l goes along the lines depicted 
in Alg0] In the first for loop, we compute N — 2 right singular 
vectors of matrices whose size is successively reduced. For 
this step, the complexity is 0{2k Ij)}- The 
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computation of vectors u. v and w has complexity 0{ (2k + 

Next, the second for loop performs N — 2 successive 
projections of the rows of matrices V n onto the vectors w. 
We need here 2 'fZiLi 2 (HjLij) operations. 


input : T € K hxl2X '" xlN : input data 

output: X £K IlXl2X "' xlN : rank-1 approximation 

v 0 <r- T(iy, 

V^V 0 - 

for n = 1 to N — 2 do 

v n 4— first right singular vector of V ; 

V n «- Unvec{v n ) € K 7 "+ lX/ '“+ 2 - r "+ 3 " /jv ; 

V^Vn, 

end 

(u, v) 4— first left and right singular vectors of V ; 

w <— v* tEl u; 

for n = N — 2 to 1 do 
X (n) (Vn- 1W) W H \ 
w = uec(X (n) ); 

end 

X(x) is the mode-1 unfolding of X. 

Algorithm 4: SeROAP algorithm 

For large dimensions and small N, the complexity of 
SeROAP is dominated by: 

r N 

O \(2k + 2)l[l j 

which can be significantly smaller than that of THOSVD. An 
example of typical executioij]] is given in the Appendix. 

3) Description and complexity of DCPD: The DCPD is 
an iterative deflation algorithm 1241 that computes the CP 
decomposition for real or complex tensors. As summarized in 
Algj5] it proceeds as follows. In the first for loop, we compute 
the rank-1 tensors X [1,1],..., X[R , 1] by successive rank-1 
approximations and subtractions. Since the rank of the tensor 
does not decrease with subtractions in general ED, a residual 
£[R, 1] is then produced. In the iterative process (repeat loop), 
a new rank-1 component is generated from the sum of the 
previous residual and Af [1,1], and a new residual £ [1,2] is 
produced with the subtraction ^[1,2] — Af [1,1]. The tensor 
y[l,2] is updated within the if-else condition. By applying 
the same procedure to the other components, we update all R 
rank-1 tensors, so that another residual £[R, 2] is generated 
in the end of the second for loop. The second loop continues 
to execute until some stopping criterion is satisfied, and all 
rank-1 components of 7” can be recovered. 

The complexity per iteration is dominated by the rank- 
1 approximation function f, which is computed R times. 
Therefore, the complexity of DCPD is 


input : T G K IlXl2X '" xlN : input data, 

R: rank parameter. 

<jy. an algorithm computing a rank-1 
approximation 

output: X r e K IlXl2X '" xlN , for r = 1, ..., R: rank-1 
components 

y[l,l]<-T; 

for r = 1 to R do 

X[r, 1] 1]); 

if r < R then 

I y[r + 1,1] t— y[r, 1] — X[r, 1]; 

else 

I £[R,l]i-y[R,l]-X[R,l]; 

end 

end 

l <- 2; 

repeat 

for r : 1 to R do 

if r > 1 then 

| y[r,l\^X[r,l-l] + £{r-l,l\-, 

else 

| y[i,i] <- x[i,i-i] + £[R,i-iy, 

end 

£[r,l]<-y[r,l\-X[r,l]\ 

end 

l <-l + 1; 

until some stopping criterion is satisfied', 
foreach r G [1,..., R] do 
| X r ^X[r,l\; 

end 

Algorithm 5: DCPD algorithm 


for the T-HSOVD algorithm, and 

N ) 

(2k + 2)Rj[lA, (4) 

for the SeROAP algorithm. 

IV. rank-1 Approximation 

This section is divided in two parts. The first one presents 
a more detailed study of the THOSVD and SeROP rank- 
1 approximation methods. For three-way tensors, we show 
that SeROAP is a better choice than THOSVD because the 
former presents a better rank-1 approximation, which ensures 
a more probable monotonic decrease of the residual ||£[i?, Z]|| 
within the DCPD algorithm l24l . A new rank-1 approximation 
algorithm is described in a second part, and is proved to 
perform better than any other rank-1 approximation method. 


#op = O ■ 


#op = 0 


N 


(2 Nk + 2)R [I Ij 


3=1 


(3) 


Matlab codes are available at http://www.gipsa-lab.grenoble-inp.fr/~pierre 
comon/TensorPackage/tensorPackage.html. 


A. THOSVD vi SeROAP 

Since we do not have at our disposal an efficient method 
to compute quickly the best rank-1 approximation of a tensor, 
we should compute a suboptimal rank-1 approximation in a 
tractable way. The THOSVD and SeROAP algorithms can 
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perform this task, as presented in Section [HI] The question that 
arises is: which algorithm performs the best? The proposition 
l4.1l below shows that SeROAP performs better than THOSVD 
for three-way tensors. For simplicity, the notation of the 
unfolding matrices does not present indices in this section. 

Proposition 4.1: Let 7” € K /lX/2X/s be a 3-order tensor. 
Let also 4> th (T) and L Sc (T) be the rank -1 approximations 
delivered by THOSVD and SeROAP algorithms, respectively. 
Then the inequality ||7~ — 0 Sc (7~) II < \\T — <f> TH (P)\\ holds. 

Proof: Let T, T™ and T £ be some mode unfolding 
of tensors r, </> TH (T) and <p Se (P), respectively. Assuming 
mode-1 unfolding for the THOSVD algorithm, we have 

||T-T ™|| 2 = ||T-Au 1 (u 3 Elu 2 ) r || 2 , 

where A, U \, u 2 , and M 3 are obtained from Alg. [3] Since A 
is the contraction of 'T on U. we plug it into the previous 
equation and we obtain, after simplifications, 

||T-T™|| 2 = ||T|| 2 HA| 2 , 

with A = u^T(ul IE u 2 ). Yet, u^T = since 

(ui.Vi, HTH 2 ) is the dominant singular triplet of matrix T. 
Hence |A | 2 = ||Tj | 2 \v ^E m^)! 2 . 

On the other hand for SeROAP, we have 

|| T - TfW 2 = || T - Tww h || 2 = 

\\T\\ 2 f + |[Tti»||l IMI 2 - 2 ||Tmi|| 2 = ||T|| 2 - w H T H Tw, 

where w is the same vector computed before the second 
for loop given in Alg. Q] for 3-order tensors. The eigenvalue 
decomposition of T H T can be expressed by 

T h T = \\T\\ 2 Vl v? + S, 


with T (x,y,a) = Ei|=i ll f i 3 - a i3 (y* E x)\\% and a = 
[«!••• a/ 3 ]. 

Plugging the optimal value of a:, 3 into Problem (O, namely 
(y* 0 x) H ti 3 , we can rewrite it as the following equivalent 
maximization problem 

[x 0 pt,y 0 P t] = argma,xz H Mz 

xy ( 6 ) 

s.t. z=y*Mx , \\y \\ 2 = 1, ||as|| 2 = 1, 

where A^ES^Vg. 

Now, we decompose M as a sum of Kronecker products. 
This can be done by reshaping M and applying the SVD 
decomposition 041 . Thus, M can be given by 

R' 

M = ^Q (r) EP (r) , 

r= 1 

with the Hermitian matrices € K /lX/l and Q^ € 

K /2X/2 . R' is the Kronecker rank of M satisfying R' < 
min(/ 2 , 1 2 ). Substituting M into Problem ©, we have: 

r! 

[xo P t,y op t\ = argma,x^ 2 (y H Q {r) *y)(x H P (r) x) 

x,y (7) 

r= 1 

s-t-Wvh = 1 , ||®||2 = l- 

Let C be the Lagrangian function given by 

£ = - ^ 2 {y H Q (r) *y){x H P (r) x)+yi{\\y\\l - l)+??2(||*||l -1), 

r =1 

where ry and 772 are the Lagrange multipliers. By computing 
the critical points, we obtain a pair of coupled eigenvalue 
problems 


where S is a semidefinite positive matrix. Hence, we have 
w h T h T-w = \\T\\\w H V\V^ w + c, 


with c > 0. To complete the proof of the proposition, we just 
need to show that \v^w \ 2 > \v^(u^ Em^)! 2 , or equivalently 
that 


l(«M>i)| > |(u^ IEm^ui)|. 


This is true, because w is by construction (cf. Alg. E0> the 
vector closest to «i among all vectors of the form a 0 b 
where a and b have unit norm. ■ 


B. Coupled-eigenvalue rank-1 approximation 

This section presents an alternating eigenvalue method for 
three-way tensors that can improve local solutions obtained 
from any other rank-1 approximation method (e.g. SeROAP 
and THOSVD algorithms). Actually, simulations have shown 
that the global solution is always attained if the initial approx¬ 
imation is close enough. 

Let ti 3 be the vectorization of slice * 3 , 1 < 13 < I 3 (we have 
chosen the third mode) of tensor 7~. The rank-1 approximation 
problem can be stated as 


[a op t,x opt ,y t ] = arg min T(x,y,a) 

^ at,x,y 

s.t. \\x\\ 2 = 1, ||y|| 2 = 1 , 


(5) 


y H A {1A) y 
_ y H ^-(h,i)y 

and 

x H B {hl) x 

. xH B {l2 ^i)X 


y H A {1Jl) y 

y H A {IlJl) y 

x H B (lj2) x 

x H B {l2j2) x 


x = Xx (8) 


y = Ay, (9) 


where A = 771 = 772 , -4( m ,n) = E PmnQ (r) *, and B (k>l) = 

r =1 
Ft' 

E Qki* p{r) ’ w ith 1 < m, n < I\ and 1 < k, l < I 2 . 

r =1 

The coupled-eigenvalue algorithm is presented in Alg. [6] 
We can initialize the algorithm by computing Xq and y Q from 
the rank-1 approximation solution obtained with SeROAP, 
THOSVD or any other rank-1 approximation method. 

The complexity per iteration of the CE algorithm is dom¬ 
inated by the construction of the matrices in the LHS of © 
and @, which is of order 0{min(/ 2 /f, I 2 ! 2 , 1 2 li)}- Suppose 
I 3 is the largest dimension. If I 3 I\ I 2 . then we can take 
advantage of the CE algorithm in terms of complexity in 
comparison with the ALS algorithm. Indeed, the complexity 
per iteration of ALS for rank-1 approximation is of order 
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input : (j>( 7 ~): rank-1 approximation 
output: ()>* {T)'. improved rank-1 approximation 
Compute tco from </>(7~) as 

x 0 4- 0(7”) (:,i 2 ,*3)/||</ > ( - 7 _ )(:,* 2 ,* 3)||2 for some i 2 ,i 3 
t 4- 0 

repeat 

Set x = x t in eigenvalue problem (9) and take y t+1 as the 
eigenvector whose eigenvalue is maximum; 

Set y = y t+l in eigenvalue problem (8) and take xt+i as 
the eigenvector whose eigenvalue is maximum; 
t 4- t + 1 

until some stopping criterion is satisfied ; 
for i 3 = 1 to I 3 do 
I «* 3 «- {U 3 ,y* t Batt); 

end 

4 - x t ® y t ® ol*\ 

Algorithm 6: CE rank-1 approximation. Above, we chose 
to start with x, but we could equivalently have started with 

y- 


(I)(.3/| 1 2 /a), which is higher than that of the CE algorithm in 
this case. Notice, however, that a properly comparison makes 
sense if the same initialization is employed in both algorithms. 

The following proposition shows that the above algorithm 
improves (in worst case the solution remains the same) any 
rank-1 approximation algorithm. 

Proposition 4.2: Let (jf'T) be a rank-1 approximation of a 
three-way tensor T. If o(fT) is the input of the CE algorithm 
and <f>*(T) the output, then the inequality \\T — <^* ('7 — )|| < 
\\T-<j>(T) || holds. 

Proof: Plugging the expression of Aj rn _ /n) into equation 
©, we obtain, after simplifications, 

a = Y,(.y H Q {r> y)i* H P {r) *), 

r =1 

which is the objective function of Problem £7}. The same result 
is obtained when the matrix is plugged into equation 

©. Now Vt > 1, let A^ and X[ v ^ be the maximal eigenvalues 
whose eigenvectors are x t and y t , respectively. 

The eigenpair (A fr+iiVt+i) obtained by solving equation © 
with x = Xt, is solution of the maximization problem 


A *+i = Q {r) *y)^. x ? p{r) x t)- 


R' 


I|W||2 = 1 


r= 1 


Also, the eigenpair (Aj + j, x t +\) obtained by setting y = y t+1 
in equation ©, is solution of the problem 


A*+i = max J2(y?+iQ {r) *y t +i)( xHp{r) x)- 

* 2 = 1 - 


R' 


r= 1 


Since 


R' 


max 
IMh = 


1 ^2(y?+iQ {r) *yt+i)( xHp{r) x) 


R' 


= ^(y?+iQ ir) *yt+i)(x?+i pir) xt+i), 


it follows in particular that 


R' 


J 2 (y?+iQ (r) *yt+i)( x ? p(r)x t) 


< f>f +1 Q«* yt+ 0(*f +1 P^Vi), 


r= 1 


which implies that A^\ < Aj*\. 

Similarly, plugging x t +i into equation ©, we can conclude 
that AS < x[% for the reason that 


R' 

r =1 

< J2(y? + 2 Q ir) *y t+ 2 )(x? + 1 P {r) x t+1 ). 

r= 1 

Hence, the sequence 

/T \ _ r \(w) \ 0 =) \(v) \( x ) \ 

i ititeN — {■ ■ •, \ 5 A t+i: A t+i> A t+2) • ■ •/ 

is monotonically non-decreasing. The same conclusion would 
be achieved if we begin by plugging x t into equation ©. 

Now, let (j>(fT) = xq ® y 0 ® c*o be a rank-1 approximation 
obtained with any other method. Assume Xq and y Q are unit 

R' 

vectors, and define Ao = {Vt) Q^*ytt){ x o p ^ x o)- By 

setting x = xq in equation (|9| in the first iteration (a similar 
operation would be possible for y 0 in equation ©), we clearly 
have Ao < A^ < A, where t max is the iteration in which 
the stopping criterion is satisfied. 

Since the optimization problems © and © are equivalent, 
a * 3 ,1 < *3 < J 3 , can be obtained by performing the scalar 
product between vectors ti 3 and yt max ^ X t m ax (which is 
equivalent to contracting tensor 7 ” on x tmax and y* ). 
Hence, the tensor 

= x tmax <g> y tmax ® a* 

is a better rank -1 approximation of 7 ~ than 4 >{T), implying 

||T-0*(T)||<||T-0(T)||. ■ 

Proposition 4.3: For any input <f>{T), the CE algorithm 
converges to a stationary point. 

Proof: In the proof of Pronosition l4.2l we have shown that 
{T t } is monotonically non-decreasing for any input o('T). 
Let p* be the maximum of the objective ©. Since the best 
rank -1 approximation problem always has a solution, then 
p * < 00 . But maxj Y(at, yt+i, a) < max^ T(x, y, a), 
which implies that {T t } t£ N is bounded above by p*. Since 
{T t } is a real non decreasing sequence bounded above, it 
converges to a limit Y*, Y * <p*. m 

V. Deflation 

In Ii24l . we proved that for a rank-i? tensor, the normalized 
residual (||£[i?, 1 } ||)jgN > 0 is a monotonically decreasing se¬ 
quence when the best rank -1 approximation is assumed within 
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Fig. 1. Visualization of the residual in an n-sphere for some iteration l of 
DCPD algorithm. 

DCPD. In this section, a thorough theoretical analysis and 
new results are presented. Based on a geometric approach, 
we sketch an analysis of the convergence of the DCDP 
algorithm, including a conjecture that it converges to an 
exact decomposition with high probability when tensors within 
‘j(R) — s^j- e j . ran k{' 7 — } < i?} are distributed according 
to absolutely continuous probability measures. 

First, let us take a closer look at the 2D geometric inter¬ 
pretation of the DCPD algorithm. Figure [T| depicts the [r,l]- 
iteration for r > 1 , so that 7 [r,l] is the angle between the 
tensors £[r — 1,1] and X[l,l — 1]. For r = 1, the residual 
£[r — 1,1] can be just replaced with £[R,l — 1] in the figure, 
and 7 [ 1 , Z] is then defined from £[R, l — 1] and X[l, l — 1], 

Before stating some theoretical results on the DCPD algo¬ 
rithm, we present a fundamental lemma related to the error in 
rank-1 approximations of tensors of the form X + £, where 
AT is a rank -1 tensor and £ any other tensor, both with entries 
in some field IK. 

Lemma 5.1: Let AT be a rank-1 tensor and (f> the best rank-1 
approximation operator. For any tensor £, 

\\X + £ - <i>{X + £)\\ <sin( 7 )||f||, 

where 7 denotes the angle between £ and X. 

Proof: Let “Px{X + £) be the orthogonal projection of 
X + £ onto span(X). Because <j>(X + £) is a best rank-1 
approximation of X + £, Vx(X + £) cannot be a strictly 
better rank-1 approximation than f(X + £). Thus, 

\\x + £ - f(X + £)\\ < \\X + £-r x (X + £)\\. 

On the other hand, X + £ — V x {X + £) _L X. Hence, we 
have \\X + £ — Vx{X + £)\\ = sin( 7 )||£|| by using basic 
trigonometry. This concludes the proof. ■ 

The following results for the DCPD algorithm stems from 
the previous lemma. 

Corollary 5.2: The inequality \\£[r, Z]|| < sin( 7 [r, l]) ||£[r— 
1, l] || holds for any 1 < r < R. 

Proof: By replacing X, £ and 7 in Lemma 15.11 with 
X[r,l — l],£[r — l,l\ and 7 [r, l] respectively, the result follows 
directly. ■ 


Corollary 5.3: For any l > 1 and q = sin('y[r', Z]), 
the inequality ||£[f?,Z]|| < ci\\£[R,l — 1]|| holds. 

Proof: By applying R — 1 times the result of Corollary 
15.21 we have 

\\£[R, l] || < (sin( 7 [R, Z]) • • • sin( 7 [ 2 , Z])) ||£[1, Z]||. 

From Lemma EU we know that 

||£[MII < sin( 7 [l,/])|| 5 [f?,Z-l]||. 

Thus, it follows that ||£[i?, Z]|| < ci\\£[R,l — 1]||. ■ 

Notice that the same result brought by Proposition (4.4) in 
[24] can be deduced from Corollary 15.31 since 0 < c; < 1 for 
every iteration l, which implies the monotonic decrease of the 
sequence {||£[i?,Z]||}z eN>0 . 

Lemma 5.4: If ||£[i?, Z]|| = ||£[f?, I — 1]||, then <7 = 1. 

Proof: Because q < 1 and ||£[iZ, Z]|| = \\£[R,l — 1]||, 
one concludes directly from Corollary 15.31 that <7 = 1. 

■ 

Lemma 15.41 shows that the DCPD algorithm might not 
improve the estimation of the rank -1 components anymore 
for l > lg > 1. And this may occur not only in the 
presence of noise. Actually, even for an almost orthogonal 
case ci « 1 , ||£[i?, Z]|| may tend to a stationary non-zero value 
as l increases. However, the DCPD algorithm converges to 
an exact decomposition if ci < C, for all Z > 1 for some 
constant C < 1. This will be subsequently detailed by means 
of a geometric approach. 

Figure |T] can also be seen as the representation of an n- 
sphere of dimension n = J 1/2 ■ • • In — 1 in K " +1 space. /3 
is half the white cone angle defined in [0,7t/2], The direction 
of the rank-1 tensor X[r,l — 1] defines the axis of the white 
cone and varies with r or l. Under a condition on f3, we can 
state an important proposition ensuring the convergence of the 
DCPD algorithm. 

Proposition 5.5: Let 7” be a tensor such that rank{7”} < 
R. An exact decomposition is recovered by the DCPD algo¬ 
rithm if and only if there exists for every (r,l) a half cone of 
angle (3 (in white in Fig. |T}, 0 < /3 < 7t/2, such that 

8 > max min y\r,l]. 

l> 1 l<r<it 

Proof: (-4=) For any iteration l > 1, take 7[ro,Z] = 
mini< r <fl 7 [r, Z]. Notice that cj < sin ( 7 [r - o, Z]) from Corol¬ 
lary [53] By hypothesis, sin ( 7 [ro, Z]) < sin(/3), which implies 
that ||£[J?,Z]|| < ci\\£[R,l — 1]|| < sin(8)\\£[R,l - 1]||. 
Because 8 is an upper bound for 7 [r, Z], Z > 1, we have 
sin(/3)||f[JZ,l]|| > ||f[JZ,2]|| =► (sin(/3 )) z_1 ||£[f?, i] I! > 
||f [i?, Z]||. Hence, when Z —> 00 , ||f [i?, Z]|| — > 0. 

(= 7 ) Let Zo be some iteration such that Zo > 1. Without 
loss of generality, assume ||£[f?, Z]|| = 0 for Z > Zo (Zo 
can be arbitrarily large). Then (||£[i?, Z]||)/gN >0 is a strictly 
monotonically decreasing sequence for 1 < Z < Zo, otherwise 
the algorithm would converge to a nonzero constant for some 
iteration smaller than Iq. Hence, for every (r, Z) we can choose 
8 , 0 < 8 < tt/ 2 such that 8 > max min 7 [r, ZL and the 

_ J>1 l<r<it L J 

proof is complete. ■ 
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As a conclusion, if for a given iteration l, all tensors £ [r, l], 
1 < r < R, fall within the gray volume depicted in Fig. Q] (the 
complementary of the white cone), then the sequence £[r,l] 
does not tend to zero. Even if this gray volume can be made 
arbitrarily small, it is not of zero measure. So the best we 
can do is to prove an almost sure convergence of the DCPD 
algorithm to the exact decomposition under some probabilistic 
conditions. 

Lemma 5.6: If tensors r are distributed within 7' R ' ac¬ 
cording to an absolutely continuous probability measure, then 
\\£[r, (]|j are absolutely continuous random variables. 

Proof: Let D = [I\ ■ ■ • In) be a specific size of n-order 
tensors, and let 7 (R) = {T G K hx ~ xl » : T C 
Because D , any tensor T within T^P is also 

distributed according to an absolutely continuous probability 
measure. Via the DCDP algorithm, each rank-1 component 
obtained in successive deflations is also in T^p. Hence, since 
the sum (subtraction) of continuous random variables does not 
affect the continuity, the residuals £[r,l) are also absolutely 
continuous random variables. Since the norm is a C° function 
in finite dimension, \\£[r, (]|| is also absolutely continuous. ■ 

For the next developments, let Zi = ||£[i?, Z]|| and define 
the following probability for some iteration L > 1: 

F l [/3\ = P [z L < sin(/3)Z L _i < . < (sin( / 3)) i_1 zPj . 

F l [P\ can be viewed as the probability that residuals fall 
within at least one of the R white cones in every iteration 
l < L. 

The following proposition ensures a reduction of Z\ by 
a factor smaller than (sin(/3 )) i_1 after L iterations with 
high probability, if a condition on the continuity of Fl[0\ is 
assumed. 

Proposition 5.7: Let L be fixed. If 3[3 o : /3o € [ 0 , 7 t/ 2 ) 
such that Fl[/3] is continuous on [/3o,7r/2], then V£ : £ G 
(0,1], 3/3 G [/3o, 7 t/ 2 ) such that F L [/3\ > 1 — s. 

Proof: Since Fffn/2) = 1 and Fl[/3\ is continuous on 
[/3 o,tt/ 2], the proof follows directly from the intermediate 
value theorem. ■ 

Although Zi, 1 < l < L, are absolutely continuous random 
variables and g m {/3) = sin m (/)) are continuous functions for 
all to > 0, the continuity of Fl[/3\ in /3 is not guaranteed 
due to the dependence of the random variables Z±,...,Zl (Zi 
depends on Z/_i). For example, for L = 2 and Z\ = 2Z-> with 
probability 1, it is easy to check that Fi [/3] is not continuous at 
/3 = 7 r/ 2 . Indeed, lim^,,.^- ^2 [/?] = 0 whereas F 2 [ 7 r/ 2 ] = 1. 

The following conjecture claims that there exists absolutely 
continuous distributions of tensors in such that the 

probability F\ [/3(Z)] tends to 1 for some function f3(l) as 
l —> 00 , and at the same time the norm of residuals tends to 0 , 
which is suitable for the convergence of the DCPD algorithm 
to an exact CP decomposition. 

Conjecture 5.8: There exists at least one absolutely con¬ 
tinuous probability measure p for tensors If within 7^ R > for 
which the following holds: 

(i). V£ : e G (0,1], V( : l > 1, 3/3 : (3 G [0, zr/2), such that 
Fi\J3] > l-£. 


(ii). V( : l > 1, 3(3 : /3 G [0, 7 t/ 2 ), such that (sin(/3 )) / ^ 1 is 
a strictly monotonically decreasing sequence converging 
to 0 . 

Subsequent computer simulations support the existence of 
a uniform probability measure /< for the entries of tensors 
within T^, such that Fp[L] « 1 for large values of L, and 
||£[i?,L]|| « 0. This reinforces our conjecture. 

VI. Computer Results 
A. Comparison between THOSVD and SeROAP 

In this section, we compare the performance of rank-1 
approximation methods SeROAP and THOSVD for different 
three-way tensor scenarios. For each case, 300 complex tensors 
whose real and imaginary parts are uniformly distributed 
in [—1,1] were generated. Figure [2] presents the difference 
between the Frobenius norms of the residuals computed as 

A0= ||T-</> th (T)|| - ||T-</> Se (T)||. 


3x4x5 tensors 3 x 4 x 20 tensors 



3 x 20 x 20 tensors 20 x 20 x 20 tensors 



0 l-:-■- o 1 -■- ; - 

0 100 200 300 0 100 200 300 


Fig. 2. THOSVD and SeROAP comparison 

We note that A cf> > 0 in all scenarios, as predicted by 
Proposition 14.11 

B. Performance of rank-1 approximations 

The tables below compare different rank-1 approximation 
methods with respect to the best rank -1 approximation, which 
was obtained from the algebraic geometric moment method 
described in l28l . Because the latter is infeasible to compute 
for high dimensions, we have focused on 2 x 2 x 2 and 3x3x3 
real tensors. 

The two iterative methods, namely ALS and CE, are initial¬ 
ized by the result obtained from SeROAP. A sample of 200 
real tensors uniformly distributed in [— 1 , 1 ] were generated 
for each of both scenarios. For comparison, we consider the 
MSE metric given by 
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where = ||T (n) - 0 m (T (n) )|| - ||T (n) - P(T (n) )\\. 

4>m( 7 _( '”' ) ) and f* (fT^) are the rank -1 approximation for 
algorithm m and the best rank -1 approximation of T' l '"\ 
respectively. 



2 x 2 x 2 tensors 

Algorithm 

MSE 

mean iteration 

THOSVD 

0.02299 

N/A 

SeROAP 

4.36155e-4 

N/A 

CE 

1.13056e-18 

6.115 

ALS 

1.03406e-13 

6.135 



3x3x3 tensors 

Algorithm 

MSE 

mean iteration 

THOSVD 

0.08386 

N/A 

SeROAP 

0.00172 

N/A 

CE 

1.70946e-4 

10.990 

ALS 

1.70949e-4 

10.435 


The results show that SeROAP is a better rank-1 approxi¬ 
mation than THOSVD as expected. For 2x2x2 tensors, 
CE attains the best rank-1 approximation. In both scenarios, 
ALS and CE converge approximately in the same number of 
iterations. 

C. Percentage of successful decompositions 
Figure [3] presents the percentage of successful decomposi¬ 
tions of rank-3 tensors for the algorithms ALS, CG with ELS, 
and DCPD. The ALS and the CG algorithms were randomly 
initialized. We have simulated DCPD with the algebraic meth¬ 
ods THOSVD and SeROAP. Noise is not considered in this 
case so that the performance is evaluated for the computation 
of an exact decomposition of 300 tensors. We consider that a 
decomposition is succeeded if the residual ||£|| < 10 -6 



3x3x3 4x4x4 5x5x5 6x6x6 7x7x7 8x8x8 

Fig. 3. Percentage of successful decomposition for rank-3 tensors. 

We note that the DCPD algorithm combined with SeROAP 
always presents a better performance than the standard ALS 
algorithm. Moreover, for higher dimensions, the percentage 
of successful decompositions is almost 100% for DCPD- 
SeROAP, which is a remarkable result, bearing in mind that 
the objective is multimodal. 

D. Convergence rate 

Figure [4] presents the performance of the algorithms in terms 
of the average of ||£[i?, 1 ] || per iteration for different values of 
the signal-to-noise (SNR) ratio for 5x5x5- rank-3 tensors. 


SNR: 40 dB SNR: 30 dB SNR: 20 dB 



Iteration 

Fig. 4. Mean of ||£[f?, l]\\ per iteration for different values of SNR. 


Again, ALS and CG algorithms were randomly initialized. 
Additive Gaussian noise is considered in our simulations. 

We note in this figure that DCPD-SeROAP converges more 
quickly than the other algorithms. For an SNR of 40 dB, 
DCPD-SeROAP attains ||£[i?, Z]|| ~ 0.01 in approximately 
100 iterations while, for the other algorithms, ||f [i?, l] || > 0.02 
for the same number of iterations. Similar results are observed 
for other SNRs. The figure also shows that performances 
become similar when the SNR is decreased. 

E. Residual vs rank 

Now, we compare the algorithms for two SNRs by varying 
the rank of 8 x 8 x 8 tensors. Again, we note the better 


8x8x8: SNR = 40 dB 




Fig. 5. Mean of ||£|| under rank variation. 


performance of DCPD-SeROAP over the competing algo¬ 
rithms. Figure[5]also shows that the combination of DCPD and 
THOSVD yields the worst results. This is expected because 
the rank-1 approximation obtained by THOSVD is not good 
enough, so that DCPD does not come at a small residual. 

VII. Conclusion 

In this paper, we presented some CP tensor decomposi¬ 
tions algorithms and provided an analysis of their compu- 
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tational complexities. Our contributions included: (i) a new 
algebraic rank-1 method, namely SeROAP, performing better 
than THOSVD for three-way tensors; (ii) an iterative rank-1 
approximation algorithm, namely CE, that refines any rank- 
1 approximation method, such as SeROAP and THOSVD, 
which converges in very few iterations; and (iii) an analysis 
of the convergence of the DCPD algorithm under a geometric 
point of view. Several computer experiments have confirmed 
the theoretical results. 

Appendix A 

Computation of rank-1 approximation using 
SeROAP 

We present an example of how SeROAP algorithm works 
for computing a rank -1 approximation of a given tensor. 

Let Tbea2x2x2x2 complex tensor whose mode-1 
unfolding is given by 

[ 1-1 0 13*0 1 ' 

1 (1) “ [o 1 -? 1 1 0 2 - 2 ? ' 

In the first for loop of SeROAP algorithm, v \ is the dominate 
right singular vector of V = Vo = T(i) given by 

—0.1717 — 0.0914*" 

0.0245 + 0.1060* 

-0.0146-0.1472?: 

-0.3189 - 0.0768? 

Vl ~ -0.6624-0.2596? ' 

-0.0914 + 0.1717? 

-0.2944 + 0.0292?: 

-0.2010 - 0.3858? 

By reshaping V\ in a 2 x 4 matrix V i we have 

_ [-0.1717-0.0914?; —0.0146 — 0.1472? 

Vl " [ 0.0245 + 0.1060? -0.3189-0.0768? 

-0.6624 - 0.2596? -0.2944 + 0.0292?' 
-0.0914 + 0.1717? -0.2010- 0.3858? ' 

In next iteration (n = 2), we compute the dominate right 
singular vector of V = V- t . Hence, 

'-0.1654 + 0.1657?;' 

0.0611 + 0.2758? 

V2 ~ -0.6190 + 0.6261? ’ 

-0.2033 + 0.2210? 

and V = V 2 — Unvec{v 2 ) is updated 

-0.1654 + 0.1657?; -0.6190 + 0.6261?' 

0.0611 + 0.2758? -0.2033 + 0.2210? ' 

The next step is to compute the vector w = v* Kit, where u 
and v are the first left and right singular vectors, respectively. 
Thus, 

0.6106 - 0.7024?! [-0.3027 - 0.0545?' 

[0.1758 - 0.3208?J ’ [-0.9481 + 0.0809?J ’ 

'-0.1466 + 0.2459?' 
-0.0357+0.1067?; 

W -0.6357+0.6166? ' 

-0.1926 + 0.2899? 


Notice that w can be viewed as a vectorization of a rank- 
1 matrix. In the end of the first iteration of the second for 
loop, w is updated so that it becomes a vectorization of a 
rank-1 three-way tensor. This is achieved by performing the 
projection of each row of matrix V\ onto w, by means of 
V-\ w, following by the multiplication with w H . Hence 

_ [-0.1900 - 0.1178? -0.0631 - 0.0611? 

A (2) " [-0.0604 + 0.0051? -0.0236 - 0.0031? 

-0.6624 - 0.1989? -0.2377 - 0.1317?' 

-0.1762 + 0.0638? -0.0730 + 0.0098?J ’ 

and 

'-0.1900- 0.1178?' 

-0.0604 + 0.0051? 

-0.0631 - 0.0611? 

-0.0236 - 0.0031? 

W ~ -0.6624- 0.1989? 

-0.1762 + 0.0638? 

-0.2377- 0.1317? 

-0.0730 + 0.0098? 

In the next iteration, we perform X\ = {V f)w H , which 
is an unfolding of a 4-order rank-1 tensor. Actually, in every 
iteration of the second for loop, X n is updated to an unfolding 
of a rank-1 tensor of order N — n + 1. Hence, 

_ [0.5373 - 0.0992? 0.1329 + 0.0653? 

A (1) “ [o.2696 - 0.1010? 0.0750 + 0.0216? 

0.1981 -0.0830? 0.0565 + 0.0140? 

0.0951 - 0.0613? 0.0306 + 0.0020? 

1.6851 + 0.1360? 0.3446 + 0.3020? 

0.8868- 0.0849? 0.2066 + 0.1249? 
0.6585- 0.0886? 0.1576 + 0.08721' 
0.3335 - 0.1067? 0.0898 + 0.0307iJ ’ 

Indeed, Xnj is the mode-1 unfolding of the rank-1 approx¬ 
imation of 'T computed by SeROAP algorithm. 
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